Import libraries

library(tidyverse)
library(ggplot2)
library(leaflet)
library(RColorBrewer)
library(scales)
library(lattice)
library(dplyr)
library(sf)
library(tidyverse)

# For pretty knitting
library(lemon)
knit_print.data.frame <- lemon_print
knit_print.tbl <- lemon_print
knit_print.summary <- lemon_print

Import and organize data sets

## Import raw Travel Time Matrix (ttm)
ttm <- read.csv('../../data/clean/ttm.csv')
df<-read_csv('../../data/score_sets/newest_long_scores.csv')
## Parsed with column specification:
## cols(
##   fromId = col_double(),
##   type = col_character(),
##   weight = col_character(),
##   nearest_n = col_character(),
##   score = col_double(),
##   pop = col_double()
## )
# convert Ids from double to factor
ttm$fromId <- as.factor(ttm$fromId)
ttm$toId <- as.factor(ttm$toId)

target_amenities <- c('gallery', 'museum', 'library or archives', 'theatre/performance and concert hall')
amenities <- read.csv('../../data/clean/vancouver_facilities_2.csv') %>% filter(type %in% target_amenities)
# clean
amenities <- amenities[,c(1,4)] # only need id and type columns
amenities$id <- as.factor(amenities$id)     # convert to factor
amenities$type <- as.factor(amenities$type) # convert to factor

# Import weight
dest_wts <- read.csv('../../data/amenity_weight/poi_index.csv')

ttm <- ttm %>% inner_join(amenities, by = c('toId' = 'id'))
head(ttm)
##        fromId toId avg_unique_time sd_unique_time    type
## 1 59150004004   10        99.76316       5.364721  museum
## 2 59150004004   15        72.48718       3.401794  museum
## 3 59150004004  157        96.69231       3.001349 gallery
## 4 59150004004 1759       106.82051       4.388213  museum
## 5 59150004004 1760        46.58974       2.642944 gallery
## 6 59150004004 1822        76.64103       3.990035  museum
### only consider the nearest amenity
nearest_1_ttm <- ttm %>%
                  group_by(fromId, type) %>%
                  summarise(avg_time = min(avg_unique_time), 
                            toId=toId[which.min(avg_unique_time)],
                            sd_time = sd_unique_time[which.min(avg_unique_time)])
## `summarise()` has grouped output by 'fromId'. You can override using the `.groups` argument.
nearest_1_ttm<-na.omit(nearest_1_ttm)
nearest_1_ttm%>%filter(type=='museum')%>%group_by(toId)%>%filter(avg_time<15)
## # A tibble: 1,984 x 5
## # Groups:   toId [85]
##    fromId      type   avg_time toId  sd_time
##    <fct>       <fct>     <dbl> <fct>   <dbl>
##  1 59150028019 museum     14.5 9570    1.71 
##  2 59150039004 museum     12   9569    0.427
##  3 59150039005 museum     13   9569    0.427
##  4 59150039006 museum     14   9569    0.427
##  5 59150040003 museum     11   9570    0.427
##  6 59150040004 museum      6   9570    0.427
##  7 59150040005 museum     10   9570    0.427
##  8 59150040006 museum     11   9569    0.427
##  9 59150040009 museum     12   9570    0.427
## 10 59150041004 museum      5   9569    0.427
## # … with 1,974 more rows

Assign level to each group

Group1: 0-15
Group2: 15-30 Group3: 30-45 Group4: 45-60

# check the histogram of avg time 
plot(hist(nearest_1_ttm$avg_time))

# create a new column useing ifelse
# between function is inclusive <=
nearest_1_ttm$time_group= with(nearest_1_ttm, ifelse(avg_time<15,"<15",
                  ifelse(between(avg_time,15,29) , "15-30",
                         ifelse(between(avg_time,30,44),"30-45","45-60"))))


nearest_1_ttm$time_group<-as.factor(nearest_1_ttm$time_group)
nearest_1_ttm
## # A tibble: 57,282 x 6
## # Groups:   fromId [14,343]
##    fromId      type                            avg_time toId  sd_time time_group
##    <fct>       <fct>                              <dbl> <fct>   <dbl> <fct>     
##  1 59150004004 gallery                             35.8 3569     1.90 30-45     
##  2 59150004004 library or archives                 35.7 9567     1.72 30-45     
##  3 59150004004 museum                              38.5 9570     1.41 30-45     
##  4 59150004004 theatre/performance and concer…     48.2 7239     2.60 45-60     
##  5 59150004005 gallery                             37.9 3569     1.60 30-45     
##  6 59150004005 library or archives                 37.5 9567     1.52 30-45     
##  7 59150004005 museum                              40.4 9570     1.44 30-45     
##  8 59150004005 theatre/performance and concer…     50.8 7239     2.21 45-60     
##  9 59150004006 gallery                             38.9 3569     1.51 30-45     
## 10 59150004006 library or archives                 38.5 9567     1.54 30-45     
## # … with 57,272 more rows
# nearest 1 export 

write_csv(nearest_1_ttm,"../../data/score_sets/nearest_1_ttm.csv")
# Import shape file data
# keep necessary columns and rows for shp
canada_dbs <- st_read("../../code/Visualizations/census2016_DBS_shp/DB_Van_CMA/DB_Van_CMA.shp", stringsAsFactors = FALSE)
## Reading layer `DB_Van_CMA' from data source `/Users/yuxuancui/Desktop/MDS/data599/w2020-data599-capstone-projects-statistics-canada-transit/code/Visualizations/census2016_DBS_shp/DB_Van_CMA/DB_Van_CMA.shp' using driver `ESRI Shapefile'
## Simple feature collection with 15197 features and 27 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 4001643 ymin: 1957237 xmax: 4068119 ymax: 2032663
## projected CRS:  PCS_Lambert_Conformal_Conic
van_dbs <- data.frame(canada_dbs[which(canada_dbs$CMANAME == "Vancouver"), ])
van_dbs$ID <- seq.int(nrow(van_dbs))
clean_van_dbs <- van_dbs[, c(1, 28, 29)]
head(clean_van_dbs)
##         DBUID                       geometry ID
## 1 59150244005 MULTIPOLYGON (((4031013 200...  1
## 2 59150244008 MULTIPOLYGON (((4031720 200...  2
## 3 59150372006 MULTIPOLYGON (((4019999 200...  3
## 4 59150372007 MULTIPOLYGON (((4019845 200...  4
## 5 59150373001 MULTIPOLYGON (((4019658 200...  5
## 6 59150373002 MULTIPOLYGON (((4019485 200...  6
# join data into a single dataframe
# convert back to sf object
van_dbs_scores <- left_join(clean_van_dbs,nearest_1_ttm , by = c('DBUID' = 'fromId'))
van_dbs_scores_sf <- st_as_sf(van_dbs_scores)
van_dbs_scores_st <- st_transform(van_dbs_scores_sf,crs = 4326)
van_dbs_scores_st
## Simple feature collection with 58136 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -123.4319 ymin: 49.00192 xmax: -122.4086 ymax: 49.57428
## geographic CRS: WGS 84
## First 10 features:
##          DBUID ID                                 type avg_time toId   sd_time
## 1  59150244005  1                              gallery 35.71795 8216 4.3222110
## 2  59150244005  1                  library or archives 13.00000 7270 0.4268328
## 3  59150244005  1                               museum 54.53846 5738 3.9725780
## 4  59150244005  1 theatre/performance and concert hall 36.07692 3080 4.3551820
## 5  59150244008  2                                 <NA>       NA <NA>        NA
## 6  59150372006  3                              gallery 19.07692 9580 1.0854200
## 7  59150372006  3                  library or archives  8.00000 5134 0.4268328
## 8  59150372006  3                               museum 18.66667 4783 1.0596260
## 9  59150372006  3 theatre/performance and concert hall 19.07692 9581 1.0854200
## 10 59150372007  4                              gallery 19.28205 9580 1.0748000
##    time_group                       geometry
## 1       30-45 MULTIPOLYGON (((-122.9625 4...
## 2         <15 MULTIPOLYGON (((-122.9625 4...
## 3       45-60 MULTIPOLYGON (((-122.9625 4...
## 4       30-45 MULTIPOLYGON (((-122.9625 4...
## 5        <NA> MULTIPOLYGON (((-122.9683 4...
## 6       15-30 MULTIPOLYGON (((-123.078 49...
## 7         <15 MULTIPOLYGON (((-123.078 49...
## 8       15-30 MULTIPOLYGON (((-123.078 49...
## 9       15-30 MULTIPOLYGON (((-123.078 49...
## 10      15-30 MULTIPOLYGON (((-123.0799 4...
polyg_subset <- van_dbs_scores_st%>%filter(type=="museum")
score_vec <- polyg_subset$time_group

# colour palette

pal_fun <- colorFactor(
  palette = c("#e30606", "#fd8d3c", "#ffe669", "#cdff5e"),
  levels = unique(polyg_subset$time_group)
)

    
  # interactive popup
  p_popup <- paste0("<strong>Accessibility:</strong>", score_vec)
    
    # Add add-ons to maps
  

  
  leaflet(data =polyg_subset) %>%
    addPolygons(
      stroke = FALSE,  # remove polygon borders
      fillColor = ~pal_fun(time_group) , # set fill colour with pallette fxn from aboc
      fillOpacity = 0.4, smoothFactor = 0.2, # aesthetics
      popup = p_popup) %>% # add message popup to each block
    addTiles() %>%
      setView(lng = -122.8, lat = 49.2, zoom = 11) %>%
    addLegend("bottomleft",  # location
              pal=pal_fun,    # palette function
              values=~score_vec,  # value to be passed to palette function
              title = "Vancouver museum Transit Accessibility") # legend title